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Path Integral Monte Carlo was used to calculate the Bose-Einstein condensate fraction at the 
surface of a helium film at T = 0.77 K, as a function of density[l]. Moving from the center of the 
slab to the surface, the condensate fraction was found to initially increase with decreasing density 
to a maximum value of 0.9, before decreasing. Long wavelength density correlations were observed 
in the static structure factor at the surface of the slab. A surface dispersion relation was calculated 
from imaginary-time density-density correlations. 

Similar calculations of the superfluid density throughout 4 He droplets doped with linear impurities 
(HCN) n are presented[2]. After deriving a local estimator for the superfluid density distribution, we 
find a decreased superfluid response in the first solvation layer. This effective normal fluid exhibits 
temperature dependence similar to that of a two-dimensional helium system. 


I. BOSE CONDENSATION AT A FREE 4 HE 
SURFACE 

It has been suggested [9] that the condensate fraction in 
the low density region near a 4 He surface is significantly 
larger than the value in bulk helium of 0.1 [7] . Varia- 
tional Monte Carlo (VMC) simulations by Lewart and 
Pandharipande[8] of small (N = 70) 4 He droplets using 
a Jastrow one-body (JOB) trial wave function give evi- 
dence for a condensate fraction which approaches unity 
as the density goes to zero in the helium surface. How- 
ever, subsequent calculations by Galli and Reatto[10] 
have shown that the condensate fraction throughout a 
helium surface computed using VMC is highly sensitive 
to the choice of trial wave function. They found that cal- 
culations performed using a shadow wave function with 
a glue term (G-SWF) have enhanced density-density cor- 
relations at long wavelengths[ll], and a maximum con- 
densate fraction of only 0.5. Significant ripplon excita- 
tions are expected to reduce the condensate fraction at 
the surface, as fewer atoms are able to occupy the zero- 
momentum state. Quantum evaporation experiments [6] 
can be interpreted as providing evidence of an enhanced 
condensate fraction. 

To avoid the problem of trial function bias and to in- 
clude finite-temperature effects, we have used Path Inte- 
gral Monte Carlo (PIMC) to calculate the density-density 
correlation functions and condensate fraction at the sur- 
face of liquid 4 He. We have also used imaginary-time cor- 
relation functions to calculate the dispersion relation of 
surface excitations in a free helium surface, and find good 
agreement with experimental thin film measurements. 


II. PIMC SIMULATION 

Our simulation system consisted of 4 He atoms, inter- 
acting pair-wise with a very accurate potential[12]. Peri- 


odic boundary conditions were used with a box size and 
initial conditions chosen to favor a double-sided film ori- 
ented perpendicular to the z-axis. To maintain a stable 
film and minimize finite size effects, we added an exter- 
nal potential determined from the long-ranged part of 
the interaction potential and the missing atoms from the 
other side of the slab, so that atoms on each of the two 
surfaces saw a potential as if they were at the surface of 
a semi-infinite slab[13]. 
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FIG. 1: PIMC density distribution, for an N = 432, T = 
0.77 K semi-infinite 4 He slab. The N — 216 results are in- 
distinguishable. The dashed line shows the effective helium 
density felt by atoms at 2 > 0. 

Most of the calculations were performed with T = 
0.77 K, with an imaginary time step of r = 1/20 K _1 . We 
performed simulations of helium slabs containing N = 
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216 and N = 432 atoms, with dimensions 24 x 24 x 17 A 
and 34 x 34 x 17 A, respectively. Lower temperature cal- 
culations were also done to determine the temperature 
dependence. 


III. THE STATIC STRUCTURE FACTOR 


rather the degree to which the G-SWF significantly over- 
estimates the effect of ripplons in a free helium surface. 
Calculations at N = 216 and T = 0.77 K agree with the 
results of Fig. 2 within statistical error, indicating that 
finite-size effects are negligible. At T = 0.38 K, we find a 
measurable decrease in the long-wavelength correlations. 
Further studies are needed to establish the temperature 
dependence. 


In order to determine the extent to which ripplons are 
present in a free helium surface, we estimated density- 
density correlations at the surface[ll] with the static 
structure factor defined as: 

S(kgi z ,z') = (p ku (z)p ku (z')) (3.1) 

where p k (z) = ^ — V e lkri 6(zi — z), N(z) is the num- 

\/N(z) t 
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FIG. 2: S(k\\,z,z) vs. fc|| throughout the surface region in an 
N = 432 4 He slab, calculated with PIMC, for T = 0.77 K. 

ber of particles in the bin at position z and k \ | is the wave 
vector parallel to the surface. This measures the correla- 
tion between density fluctuations at vertical positions z 
and z! . 

Fig. 2 shows S(k\\,z, z) curves as a function of density 
(z). Each curve is the average of both sides of the slab 
of eight identical simulations (a total of 7 x 10 5 Monte 
Carlo passes). At densities of p{z) = 0.015 A -3 and 
below there is a small enhancement of long-wavelength 
density-density correlations at k\\ = 2i i/L = 0.18 A -1 , 
which is evidence for ripplons. However, the curves are 
closer to the VMC calculations of Galli and Reatto which 
used a JOB trial wave function than those which used 
the G-SWF form. This does not imply that the JOB 
trial wave function is well-suited to representing an in- 
homogeneous helium system such as a helium slab, but 


IV. THE EXCITATION SPECTRUM 


The excitation spectrum can be estimated with path 
integrals using imaginary-time correlation functions [14]. 
The dynamic structure factor is related to the imaginary- 
time density-density correlation function by: 

OO 

F(k,t) = f duj e~ tw S(k, w) (4.1) 

— OO 

= £ (Pk(i)Pk(O)) . (4.2) 

To select out the excitations at the free surface, we want 
to calculate the imaginary-time correlation function of 
propagating surface modes. W. F. Saam[15] proposed 
that the lowest quantized hydrodynamic mode (capillary 
wave) at a free helium surface will have the form 

<t>h o(r, t) = (j>ho{z) e* k r n e ~ lu>kot (4.3) 

where 

<MA) oc e~^ z (4.4) 

and the decay constant n{k) is defined as 

n{h) = —bk + {k 2 + b 2 ) 1 / 2 (4.5) 

h = (4.6) 

where do is the zero-temperature surface tension, po is 
the bulk density, and s is the zero-temperature sound 
velocity. To calculate the dispersion relations for excita- 
tions of this form, we use pj. in Eq. (4.2), defined as 


Pk = e ‘ k "' r<11 <t>k{z%) (4-7) 

i 

where 

0 fc (2) = e -«( fe )G-Zc) (4.8) 


and n(k) is defined by Eq. (4.5). We defined z c as the 
point in the surface at which the average density was 
equal to 10% of the bulk density. 

Extracting the dynamic structure factor by inverting 
Eq. (4.1) is ill-conditioned in the presence of statistical 
noise. It has been shown[16] that a maximum entropy 
method greatly increases the numerical stability. In the 
maximum entropy method, S(k,u>) is calculated by min- 
imizing the function 


F{S, a) 


g-(l/2 )Q(S) e a( (S) 

^ x 

Z Q Z C 


(4.9) 
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where Q (S) is the ’’likelihood” of the PIMC data given 
an S and £(S) is the entropy of a given S(k,u>) defined 
with respect to some default model, with a an adjustable 
parameter (also optimized). 


tance Zq from the center of the slab is given by 

« k n(zo) = j dr -|| e _ * k|l ' r M n(r M ; ^ 0 ) (5.1) 
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FIG. 3: Dispersion relation of free surface excitations. The re- 
sults calculated from PIMC and maximum-entropy inversion 
(filled diamonds) for an N = 216, T = 0.77 K 4 He slab are 
compared with the experimental thin film data (open circles) 
of Lauter et al. and the DFT results (solid line) of Lastri et 
al. Also shown is the dispersion relation of bulk 4 He at SVP 
(dashed line). 


The dispersion energy for a given value of k can be 
determined from the position of the maximum value of 
S(k,ui). Boninsegni and Ceperley[17] found that the po- 
sition of the main S(k,u>) peaks for liquid helium agree 
quite well with experiment, despite significant broaden- 
ing of the excitation spectrum caused by the maximum 
entropy procedure. The dispersion energy of the surface 
excitations as estimated using this procedure are shown 
in Fig. 3. The two lowest data points, at k < 0.5 A, had 
significant fitting error and are only qualitatively reliable. 
Otherwise, we see excellent agreement with the experi- 
mental thin film data of [3], including the curvature of 
the ripplon branch toward the roton minimum, proposed 
as evidence for roton-ripplon hybridization[18, 19]. 


V. THE BOSE-EINSTEIN CONDENSATE 
FRACTION 


where the off-diagonal single particle density matrix is: 


n( r; z) 


1 

p{z)Z 


dr i • • • dr N p( rq, r 2 , • • • r N , ri+r, r 2 , • • • r N - (3). 


•5.2! 

and p is the many-body density matrix, Z = Tr(p) is the 
partition function, and rj| = |r || — rj||. This function can 
be calculated from PIMC [20] by performing simulations 
with a single open path. We fix the endpoints of the 
open path at z = zq, and calculate the distribution of 
end-to-end distance n(rir,Zo). The condensate fraction 
at a given point in the surface is: 


n 0 {z 0 ) 


n(r || -> oo; z Q ) 
n(r\\ -» 0; zq) ' 


(5.3) 


Using PIMC, we calculated n(r\\,z) throughout the 
slab. Nonlinear least-squares fitting was used in the re- 
gion of t*| | < 1.5 A to get an estimate of n{r\\ = 0,2), 
and n(r|| — > 00 , 2 ) was calculated by averaging over the 
region at large r|| where n(r\\,z) is flat. At the low- 
est densities, it is not clear whether the n(r)’s have 
reached their asymptotic limit within the finite simula- 
tion box. Thus, for those densities, ( for N = 216 below 
p{z) = 0.001 A" 3 ) our results are upper bounds to the 
condensate fraction. 

The condensate fraction uq{z) is plotted as a function 
of average density p(z) in Fig. 5, for both N = 216 and 
N = 432 helium slabs. As one moves through the surface, 
the condensate fraction initially increases with decreas- 
ing density, due to the decreased zero-point motion from 
helium-helium interactions, reaching a maximum value 
of 0.93(3) at p{z) = 0.002 A -3 . As the average den- 
sity decreases below this point, the condensate fraction 
begins to decrease, further evidence for correlated den- 
sity fluctuations due to ripplons at the surface. This is 
in qualitative agreement with the G-SWF VMC calcu- 
lations of Galli and Reatto. However, the G-SWF trial 
wave function significantly overestimates the degree to 
which ripplons are present in the surface. 

Our results support the model of a free helium surface 
with ripplons, in which the condensate fraction reaches 
a maximum at an intermediate density in the liquid- 
vacuum interface, before decreasing at lower densities. 
Experimental probes of the surface will indeed see an en- 
hanced condensate fraction as proposed by Griffin and 
Stringari [9]. 


We define the condensate fraction in the slab geometry 
by the fraction of atoms at a given value of 2 having pre- 
cisely k\\ =0. (Because [k\\,z] = 0, we can measure the 
momentum parallel to the surface simultaneously with 
the z-position.) The momentum distribution at a dis- 


VI. LOCAL SUPERFLUIDITY DENSITY 
ESTIMATOR 
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In this section, a microscopic Path Integral Monte 
Carlo (PIMC) estimator of the local contribution to the 
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FIG. 4: n(r\\,z) vs r||, throughout the surface region of N = 
216 (top) and N = 432 (bottom) 4 He slabs at T = 0.77 K, 
calculated with PIMC. The slab dimensions were 24x24x 17 A 
and 34 x 34 x 17 A respectively. Density labels correspond 
to z = 0,8,9,10,10.5,11,12,13 A for IV = 216, and z = 
0,9, 11,12, 13 A for IV = 432. 


total superfluid response is presented. The superfluid 
density can be defined in terms of the response of the sys- 
tem to an imposed rotation. In imaginary-time path in- 
tegrals, it is manifested by particle exchange over length 
scales equal to the system size. Although superfluid re- 
sponse (like conductivity) is not itself a local property, it 
is possible to calculate a local contribution to the total 
response. In PIMC, the total superfluid response along 
the axis of rotation n is proportional to the square of 



FIG. 5: Condensate fraction vs. density throughout the sur- 
face region of a T = 0.77 K 4 He slab. Also shown are the 
VMC calculations of Galli and Reatto. 


total projected area of the imaginary-time paths [24]: 


P 


n 


2m ( A l) 

/3A Ic ’ 


(6.1) 


where A = a /2m, (3 = 1 /ksT and I c is the classical 
moment of inertia. To define a local superfluid density 
we write: 


Ps( r )U = (f dr ' An(r)An(r')) 

f^(A*(r)Aa), (6.2) 

where A(r), related to the local angular momentum op- 
erator, is defined as 

A(r) = i T> it » x r i)fc+ i) S(r - r i)fc ) (6.3) 

i,k 

and is the n-component of the total area of all the 
particles. Since Ar integrates to the total area, the lo- 
cal superfluid response exactly integrates to the total re- 
sponse. 

Two types of contributions to the local superfluid den- 
sity can be distinguished based on the connectivity of 
the instantaneous paths: contributions of particles on 
the same chain, which on the average must be positive, 
and contributions of particles on different exchange cy- 
cles. By reversing the order of one exchange cycle the 
contribution from different cycles will change sign; if the 
cycles are spatially separated, the magnitude of the con- 
tribution will be unaffected, so that their contribution is 
much smaller. They will however increase the statistical 
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In order to test our estimator, PIMC simulations were 
performed on N = 128 and N = 500 4 He droplets 
doped with (HCN) 3 isomers. Nauta and Miller [32] found 
that HCN molecules in helium droplets self-assemble into 
linear chains spaced roughly 4.4 A apart. The HCN 
molecules in our simulations were fixed along the z-axis 
with this spacing. An imaginary time-step of 1/20 K -1 
was used. With cylindrical symmetry and precise exper- 
imental data over a range of isomer lengths, this system 
is well-suited for studying superfluidity at a molecular 
interface. 

The number density of a doped helium droplet at 
T = 0.38 K is shown in Fig. 6. The two-dimensional 
areal density of the first solvation layer was 0.12 A -2 , 
which is still in the liquid phase for thin 4 He films at 
these temperatures [25]. Even if the density were large 
enough to be solid in a 2D film, the curved geometry of 
the film around the HCN molecule could frustrate crys- 
talline order, keeping the first solvation layer in a dense 
liquid state. 


FIG. 6: Average number density of an N = 500 4 He-(HCN) 3 
droplet at T = 0.38 K. The grey scale saturates at p = 
0.03 A -3 (black). The dashed lines define the cylindrically- 
symmetric region used for averaging superfluid density. 

noise of the superfluid density. Numerically we find that 
the same cycle contribution accounts for roughly 80% of 
the total superfluid density in the systems presented in 
this paper. 

The local superfluid density estimator used by Kwon 
and Whaley [23] defined the effective normal fluid induced 
by the anisotropic molecule-helium interaction potential 
in terms of the average number density distribution of 
paths in exchange cycles of fewer than six atoms. Al- 
though qualitatively interesting, this estimator is not the 
superfluid response to an imposed rotation, as there is 
no direct relation between the number of atoms in a per- 
mutation cycle and its area. 



VII. SUPERFLUIDITY OF HELIUM DROPLETS 
WITH MOLECULAR IMPURITIES 

Molecules confined in helium nanodroplets have been 
shown to exhibit excitation spectra with clearly re- 
solved rotational fine structure consistent with that of 
a free rotor, though with an increased moment of inertia. 
Grebenev et aZ. [21] have shown that only a few layers of 
4 He coating the molecule are required to decouple the im- 
purity from the droplet and achieve free rotation. Calle- 
gari et al. [22] have suggested that the increased moment 
of inertia is due to the hydrodynamic response of the 
impurity rotating through the anisotropic helium den- 
sity immediately surrounding the molecule. Kwon and 
Whaley [23] have proposed a model in which a micro- 
scopic normal fluid is induced in the first solvation layer 
by the anisotropy of the molecule-helium interaction. 


FIG. 7: Superfluid and number density distributions, aver- 
aged over the cylindrically-symmetric region of an N = 500 
4 He-(HCN) 3 droplet at T = 0.38 K. 

We define the change in the moment of inertia due to 
the helium as the contribution due to the normal fraction 
in the first layer rotating rigidly with the impurity: 

A 1= f drm 4 r 2 ± (p(r) - p s (r)) (7.1) 
•MS) 

where v(£) is the volume of helium a distance £ away 
from the surface of the molecule, rj_ is the radial distance 
from the axis of rotation in cylindrical coordinates. Both 
the estimated statistical error and the uncertainty due 
to the cut-off were on the order of 10%. The moment of 
inertia due to the effective normal fluid in the first sol- 
vation layer calculated using our PIMC local superfluid 
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density estimator is in reasonable agreement with the ex- 
perimental value of A / = 1240 amu A 2 [31], within error 
bars. We also find that A / is effectively independent of 
temperature below 

T = 1.0 K. This is because the dominant contribution 
to the moment of inertia comes from the induced normal 
fluid at the ends of the isomer, which is almost completely 
due to anisotropy in the molecular potential. 



V A ) 
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FIG. 8: Superfluid density fraction distributions, averaged 
over the cylindrically-symmetric region of an IV = 128 4 He- 
(HCN )3 droplet at T = 0.19, 0.38, 0.77, and 1.25 K. The top 
graph shows superfluid response about the radial axis, the 
bottom graph shows superfluid response about the molecule 
axis. 

To quantify the reduced superfluid response in the first 
layer, we averaged over the cylindrically-symmetric re- 


gion of the 4 He-(HCN) 3 droplet, defined as the region 
between z = —3.5 A and z = 3.5 A (see Fig. 6). The 
results of this averaging are shown in Fig. 7. Taking 
the ratio of the superfluid density to the number den- 
sity clearly shows a reduced superfluid response in the 
first layer for rotation about both the radial axis and the 
molecular axis. At zero temperature, there can be no 
“normal” density for rotation about the molecular axis 
because of the cylindrical symmetry, so this reduction 
must be due to thermal excitations. 

Shown in Fig. 8 is the temperature dependence of the 
superfluid density in the first solvation layer, determined 
for N = 128 4 He-(HCN) 3 droplets at T = 0.19, 0.38, 
0.77, and 1.25 K. At first glance, this result appears to 
contradict known properties of liquid helium; the super- 
fluid density of bulk helium at 0.75 K is 1.000(1). The- 
oretical studies of pure N = 128 helium clusters show a 
superfluid transition roughly in agreement with the bulk 
lambda transition [26]. Pure droplets like those produced 
for use in scattering experiments, with several thousand 
atoms at T = 0.38 K, should have a superfluid fraction 
very close to unity. However, the helium in the first 
layer coating the impurity molecule more closely resem- 
bles a two-dimensional system than a three-dimensional 
system, because the motion of the helium atoms is re- 
stricted by the He-HCN potential to be primarily on the 
cylindrical surface circumscribing the (HCN) 3 impurity. 
Two-dimensional helium films have been shown to have a 
Kosterlitz-Thouless type of superfluid transition at tem- 
peratures significantly lower than T\ [27, 28] . The reduc- 
tion in the transition temperature is due to the reduced 
dimensionality, increasing the phase space for long wave- 
length fluctuations, and the lowering of the “roton” gap. 
A similar temperature-dependent reduction in the super- 
fluid response of the first layer of helium surrounding the 
ends of the isomer was not observed. 

To extract the average superfluid density fraction in 
the first layer from the distributions plotted in Fig. 8, we 
integrate over the superfluid density and number density 
in the first layer and take the ratio: 

r i 

/ dr cy i p s {r C yi) 


f dr cy i p(r cy i) 
o 

where rq is the position of the density minimum between 
the first and second solvation layers. The average su- 
pcrfluicl density in the first solvation layer as a function 
of temperature is shown in Fig. 10. The transition is 
very broad due to the small number of atoms (~ 30) 
in the first layer. It shows the onset of superfluidity at 
roughly 1 K. The superfluid density about the molecule 
axis (z-axis) goes to unity as the temperature goes to 
zero, as the density is symmetric about this axis. The 
normal fluid response from rotating about the radial axis 
(r-axis) combines both thermal effects and effective nor- 
mal fluid induced by rotating through the anisotropic 
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helium density, and because of this, does not go to zero 
at T = 0 K[22]. 

We have determined that only a small fraction of the 
particles in the first layer are localized (not permuting) 
at T=0.38 K and below, as shown in Fig. 9. Though 
many of the permutations are between atoms within the 
same layer, the first layer is not cut off from the rest of 
the fluid. Below 1 K, the majority of the atoms in the 
droplet are part of exchange cycles with atoms in both 
the first layer and the rest of the droplet. However, in 
terms of excitations, the first layer is well represented as 
a 2D superfluid. 



FIG. 9: Exchange behavior of the first layer (defined as 2.5 A 
or less from the surface of the (HCN )3 isomer) as a function of 
temperature. A permutation cycle was considered to be part 
of the first layer if more than M/2 consecutive time slices were 
contained within it, where M is the number of imaginary time 
slices per atom. 

To compare the effects of density on the thermally- 
induced normal fluid in the first layer, we calculated the 
superfluid density distribution for N = 128 4 He-(HCN) 3 
droplets at T=0.19, 0.38, 0.77, and 1.25 K, with a 4 He- 
HCN interaction reduced by a factor of two in order to 
reduce the density in the first solvation layer. The den- 
sity in the first layer decreased by ~ 15%, correspond- 
ing to an average coverage of 0.10 A -2 . The difference 
in coverage caused a 20% reduction in the superfluid re- 
sponse about the radial axis at T=0.19 K and T=0.38 K. 
At higher temperatures, the superfluid response was un- 
changed within the estimated error bars. The super- 
fluid response about the molecule axis (z-axis) was un- 
changed within error bars. This is further evidence that 
the normal response in the first layer is due to both the 
anisotropy of the molecular potential and thermal exci- 
tations. 


Using PIMC and a local superfluid density estima- 
tor, we find that the first solvation layer surrounding 
an (HCN) 3 isomer has significant thermal excitations at 
temperatures as low as T=0.19 K, with a superfluid tran- 
sition temperature similar to that of a two-dimensional 
system. It remains to be seen how strongly these thermal 
excitations are present in droplets with different dopants; 
this is the first observation of thermal excitations in a he- 
lium droplet at such low temperatures. Although we find 
this effect only makes a small contribution to the moment 
of inertia of this system, it may explain the anomalous 
Q branch observed in helium droplets doped with long 
linear impurities. The role of thermal excitations could 
also be tested by varying the temperature of the droplet, 
by adding 3 He atoms. 



FIG. 10: Superfluid density fraction of the cylindrically- 

symmetric portion of the first layer of an N = 128 4 Hc- 
(HCN) 3 droplet vs. temperature, calculated using Eq. (7.2). 
The first layer in these systems has an estimated coverage of 
a — 0.12 A -2 . Superfluid response about both the molecule 
axis (solid squares) and radial axis (solid circles) are pre- 
sented. Also shown are the PIMC results of Gordillo et al 
(open diamonds), for a 2D 4 He film with a = 0.051 A -2 . For 
reference, the superfluid transition in bulk 4 He is included 
(dashed line). 

The authors would like to acknowledge K. Lehmann 
and B. Whaley for useful discussions. The computations 
were carried out at NCSA and the IBM cluster at the 
Materials Computation Center, and was supported by 
the NASA Microgravity Research Division, Fundamental 
Physics Program. This work was also performed under 
the auspices of the U.S. Department of Energy by Uni- 
versity of California Lawrence Livermore National Labo- 
ratory under contract No. W-7405-Eng-48. 


84 





[1] E. W. Draeger and D. M. Ceperley, Phys. Rev. Letts. 89, 
015301 (2002). 

[2] E. W. Draeger and D. M. Ceperley, Phys. Rev. Lett. 90, 
065301 (2003). 

[3] H. J. Lauter, H. Godfrin, V. L. P. Frank, and P. Leiderer, 
Phys. Rev. Lett. 68, 2484 (1992). 

[4] D. O. Edwards and W. F. Saam, in Progress in Low Tem- 
perature Physics VII A, edited by D. F. Brewer (Nortli- 
Holland, Amsterdam), p. 283, 1978. 

[5] A. Lastri, F. Dalfovo, L. Pitaevskii, and S. Stringari, J. 
Low Temp. Phys. 98, 227 (1995). 

[6] A. F. G. Wyatt, Nature 391, 56 (1998). 

[7] O. Penrose and L. Onsager, Phys. Rev. 104, 576 (1956). 

[8] D. S. Lewart, V. R. Pandharipande, and S. C. Pieper, 
Phys. Rev. B 37, 4950 (1988). 

[9] A. Griffin and S. Stringari, Phys. Rev. Letts. 76, 259 
(1996). 

[10] D. E. Galli and L. Reatto, J. Low Temp. Phys. 113 , 223 
(1998). 

[11] D. E. Galli and L. Reatto, J. Phys. Condens. Matter 12, 
6009 (2000). 

[12] R. A. Aziz, A. R. Janzen, M. R. Moldover, Phys. Rev. 
Letts. 74, 1586 (1995). 

[13] M. Wagner and D. M. Ceperley, J. Low Temp. Phys. 94, 
161 (1994). 

[14] D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995). 

[15] W. F. Saam, Phys. Rev. B 12, 163 (1975). 

[16] J. E. Gubernatis, M. Jarrell, R. N. Silver, and D. S. Sivia, 
Phys. Rev. B 44, 6011 (1991). 

[17] M. Boninsegni and D. M. Ceperley, J. Low Temp. Phys. 
104 , 339 (1996). 

[181 E. Krotscheck, S. Stringari, and J. Treiner, Phys. Rev. B 
35 , 4754 (1987). 

[19] L. Pitaevskii and S. Stringari, Phys. Rev. B 45, 13133 
(1992). 

[20] D. M. Ceperley and E. L. Pollack, Can. J. Phys. 65, 1416 


(1987). 

[21] Grebenev, S., J. P. Toennies and A. F. Vilesov, Science 
279 , 2083 (1998). 

[22] Callegari, C., A. Conjusteau, I. Reinhard, K. K. 
Lehmann, G. Scoles and F. Dalfovo, Phys. Rev. Lett. 
83 , 5058 (1999). 

[23] Kwon, Y. and K. B. Whaley, Phys. Rev. Lett. 83, 4108 
(1999). 

[24] Pollack, E. L. and D. M. Ceperley, Phys. Rev. B 36, 8343 
(1987). 

[25] Bishop, D. J., J. E. Berthold, J. M. Parpia and J. D. 
Reppy, Phys. Rev. B 24 , 5047 (1981). 

[26] Sindzingre, P., M. L. Klein and D. M. Ceperley, Phys. 
Rev. Lett. 63 , 1601 (1989). 

[27] Ceperley, D. M., and Pollock, E. L., Phys. Rev. B 39, 
2084 (1989) 

[28] Nyeki, J., R. Ray, V. Maidanov, M. Siqueira, B. Cowan 
and J. Saunders, J. Low Temp. Phys. 101, 279 (1995). 

[29] Gordillo, M. C. and D. M. Ceperley, Phys. Rev. B 58 , 
6447 (1998). 

[30] Kim, K. and W. F. Saam, Phys. Rev. B 48, 13735 (1993). 

[31] Nauta, K., private communication (2001). 

[32] Nauta, K. and R. E. Miller, Science 283, 1895 (1999). 

[33] Goyal, S., D. L. Schutt and G. Scoles, Phys. Rev. Lett. 
69 ,' 933 (1992). 

[34] Pi, M., R. Mayol and M. Barranco, Phys. Rev. Lett. 82 , 
3093 (1999). 

[35] Lehmann, K. K., Mol. Phys. 97 , 645 (1998). 

[36] Toennies, J. P. and A. F. Vilesov, Annu. Rev. Phys. 
Chem. 49, 1 (1998). 

[37] Atkins, K. M. and J. M. Hutson, J. Chem. Phys. 105, 
440 (1996). 

[38] Kwon, Y., P. Huang, M. V. Patel, D. Blume and K. B. 
Whaley, J. Chem. Phys. 113 , 6469 (2000). 


85 



